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The Telecommunications Division has built a stability analyzer for testing Deep 
Space Network installations during flight radio science experiments. The low- 
equency part of the analyzer operates by digitizing sine wave signals with band- 
widths between 80 Hz and 45 kHz. Processed outputs include spectra of signal 
phase amplitude, and differential phase; time series of the same quantities; and 
Allan deviation of phase and differential phase. This article documents the digital 
signal-processing methods programmed into the analyzer. 


I. Introduction 

a nH The , F fi Cei ! tly d 7° lOPe i radi ° SCienCe Stability anal y zer ( RSA ) 18 an instrument for real-time testing 
and certification of Deep Space Network (DSN) equipment to be used during gravity wave and planetary 

Td 1 ' TW ° Ta T P 7T ^ ^ t6Sted: (1) thG radi ° SdenCe ° pen - lo °P reiver 

s^tem fFTS ^ distribution network * the DSN frequency and timing 

sL wave signals The ^ tw0 . SOUrces are downconverted to low-frequency band-limited 

to video mfvi T a . ge ° f the °P en 'i°op receiver, called radio science intermediate frequency 
to video (RIV), produces sine wave signals with frequencies and bandwidths ranging from 150 Hz in 
an 82-Hz band to 275 kHz in a 45-kHz band; these depend on the choice of RIV fiber RIV signal 
are processed directly by the low-frequency RSA circuitry. Pairs of 100-MHz FTS signals are processed 
portion of the RSA called the 100-MHz interface assembly (100 MHz IA), which resides near the 

ThG 100 ^ mbCeS tHe tW ° Signals at 10 GHz -d downcon^rts the mi^r 

btCue^sTcilcuter ‘ b “ dWidth ' " h ‘ Ch ‘ S OV " “ fiber -° PtiC *° tl " 

i„foT!lto ft X n Th CirCUlt T 1,85 t W ° mMb “ dS “" verti "g * ^d-ltoited sine wave signal to digital 
information. First the signal can be sampled with a 16-bit analog-to-digital (A-D) converter clocked 

,7 . a synthesizer. In this mode, two signal channels can be accommodated with the aim of extracting 
heir differential phase. The maximum total data rate is about 230 kilosamples per second. Second if 

* !. "" er , TTl " f SP rS Xim ™ ly 01 Hz> " C “ ** mixed With the ° ut P« ° f onothet 

synthesizer set to this frequency minus 1 Hz. The 1-Hz mixer output is filtered and hard limited by a 

zero-crossing detector, and the up-crossing times of the resulting sequence of pulses are captured by a 
time-interval counter according to the “picket fence” method [4] captured by a 

The principal aim of processing the A-D data is to reduce their bandwidth by a user-selected factor 
and to extract the amplitude and phase modulations that constitute the sidebands of the sine wave signal.’ 
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The phase of two channels can be combined into differential phase. Three output types can be generated, 
spectrum of the signal and its modulations, time series of the modulations, and Allan deviation of P • 
As described below the digital signal processing operates in three alternate modes, called full band, 
medium band and narrow band. The choice among these depends on the desired bandwidth reductio 
factor. The 1-Hz zerocrossing data are processed in the same way as sequences of phase residua 
produced by narrow-band processing. 

The digital signal processing (DSP) methods are designed to take advantage of the architecture of 
. a^Sot vector^ processor based on the 40-MHz Intel I860. Most of the heavy hftmf ns d«e by 
manufacturer-supplied vector library routines, which include fast Fourier transform ( ) a 

pulse response (FIR) filtering routines. Throughputs of approximately 25 million floating-point operatio 

per second were achieved. 

The remainder of this article explains the DSP methods in some detail. 


II. Signal Properties 

A. Radio Frequencies 

In any test setup, there are two radio fluencies of jnterest^Let IndMan 

primary comparative mmng takes ^ b ^ ^ MHz (X -band). For an FTS test, 

deviation^ For a RIV test f J e of the 100 -kHz output of the 100 MHz IA is 

fppro^imately^ S/Snes the difference between the phases of the two 100-MHz inputs. Phase results are 
scaled by /ref //mix- 

B. Analog Sine Wave Signal 

The downconverted signal is assumed to lie in an analog frequency band with the center at / 0 fst and 
The downcon g ri V filter or the 100-MHz IA. The frequency / of ,t can 

^ posS or nSktle si 2 discussion of polarity below. Somewhere in this band - the earner. 
Fveent in full band processing, it is assumed that the signal consists of a carrier with weak sidebands, 
fhe total carrier-to-noise ratio should be at least about 30 dB. (This instrument is a stability analyzer, 

not a receiver.) 

C. Digitized Sine Wave Signal 

The analog signal is sampled by a 16-bit A-D converter at the sample rate /., which l has to be ^" 
fhat the analog frequency band is aliased into the Nyquist band (0,/ s /2) or (-/ s /2,0). In this y, 
Lhletnds J thfe-ie) ™ .—-I Each RIV H. is Signed fo, a certain In any case, an 

acceptable f, can be obtained from the formulas 


intf^i-0.5 

V vVvid 


4 1/ofstl 
2m + 1 


rrssss szzszsz: 


a = 0.944272, 


4 |/ofstl 
2m + a 
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The number a is related to the golden ratio (\/5 - l) /2. 

D. Polarity 

In the radio science receiver, the 2.3-GHz or 8.4-GHz signal is downconverted and filtered three times 
until the carrier is at 10 MHz -f / 0 f s t> where / 0 f st can be positive or negative. At this point, the spectrum 
or phase polarity of the signal is positive, i.e., the same as the radio frequency (RF) signal. The fourth 
downconversion by the 10-MHz local oscillator and subsequent filtering, therefore, yield a signal whose 
polarity equals the sign of / 0 f s t- Moreover, the sampling can flip the polarity again. To make better sense 
of this, it is good to think about the two-sided representation of the signal. One side of the signal has 
the right polarity (positive), and the other side has the wrong polarity. If we let 


^base — nint 



^pol — sign (/ofst ^bas e/s) 


where nint (x) is the nearest integer to x , then s po | is the polarity of the digitized signal, the side of 
the analog signal with the right polarity lies between n base / 5 and (n base + s po i/2) f s , and the side of the 
digitized signal with the right polarity lies between 0 and s po \f s /2. The user has the responsibility of 
entering f Q f st with the correct sign. 


III. Full-Band Processing 

This mode allows the user to see a snapshot of the signal in the time and frequency domains before 
proceeding to a closer view. The user selects an FFT size N (2048 or 4096). A frame of A-D data 
x[0], * ■ ■ , x[N - 1] is collected. These can be plotted against elapsed time in the frame, after scaling them 
back to volts at the A-D input (10 V — 32,768). A spectral estimate of the frame is computed by scaling 
the frame so that £^x[n ] 2 = 1 and calculating 


S x [k] = 


2 


^ x[n]uo[n; N, 5] exp (- i2nnk/N ) , 

n=0 


* = 0, • ■ • , N/2 


( 1 ) 


where u 0 [n; N, 5] is the Oth-order, iV-point “trig prolate” data taper [5] with bandwidth parameter w = 5 
(Appendix B), scaled so that J^uo[n] 2 = AT. The sidelobes of this taper (f^os in Fig. 1) are low enough 
so that no leakage from the carrier should be visible in the sidebands. The array 10 log i 0 5 x [fe] (labeled 
dBc/Hz) is plotted against the frequency array 


f[k\ — f s (n base -b s po \k/ N ) , k — 0, * ■ ■ , N /2 

which shows the side of the signal with the correct polarity. The user chooses how many of these frame 
spectra are averaged into a run spectrum. The frames do not have to be adjacent; it is all right to lose 
data while processing the previous frame. 

The resolution bandwidth of the spectral estimate, given by 




f s N 

(E«oN ) 2 


273 




Fig. 1. Spectral windows: full band ftos, 
medium band ^04, and narrow band O4. 

has two purposes: (1) It gives the user a rough idea of the resolution of the spectral plot, and (2) it allows 
the user to estimate the power of a bright line (narrower than W nb ) in dBc by adding 10 log xo^nb to the 
dBc/Hz reading at the peak of the line. 

Because the main purpose of this function is a check on what sort of signal is actually in the Nyquist 
band, it might be preferable to scale the spectrum to dBm/Hz or dBV 2 /Hz instead of scaling the frame 
to power 1 and claiming that we are seeing dBc/Hz. Then, for example, if no signal were present, the 
display would show the correct spectral density level of the noise. 


IV. Medium-Band Processing 

In this mode of processing, we assume that the sampled signal consists of a carrier with weak sidebands. 
The purpose of the processing is to reduce the bandwidth of the signal by a modest amount (up to 128 
with current parameters), remove the carrier, and measure properties of the sidebands. 


A. z-Frame Production 

The user having selected an FFT size Nftt and a decimation factor r, both powers of 2, define the frame 
size N x( = rNff t . In order to limit memory usage, the frame is divided into n bf adjacent batches of size 
N xb , a divisor of N x f that is not more than some maximum batch size (currently 8192). One batch at a 
time is processed. We use the first batch to measure the carrier frequency by a simple vector computation 
called “Pony, Part 1” (Appendix A). Let 6 be the measured frequency in radians per sample, the sign of 
6 being s poU and let u = exp(-io). Let *[n], n = 0, • • • , N xf - 1, be the A-D x-frame. A complex z-frame 
z[m] of size N z f < A r ff t is computed by 

zi[n] = :r[n]tt~ n , n = 0, • • • , AT xf - 1 (2) 
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n = 0, • • • , N z[ - 1 


( 3 ) 


n h — 1 

z l n ] = 'y ] [rn + A:], 


Ar— 0 


here Ar is a lowpass FIR filter designed for decimation by r (Appendix B). Its length n h is assumed to 
e a multiple of r (currently 16r), and it follows that we can take N zf = N m - n h /r + 1. The ripples of 
the frequency response of h r above the decimated Nyquist frequency (Fig. 2) are low enough so that the 
aliased image of the wrong side of the carrier at -6 barely appears above the 16-bit quantization noise 
in a spectrum output with simulated data. 

The computation in Eqs. (2) and (3) is carried out batch by batch, the z-frame being built up in n bf 
steps by an overlap-add operation. The result is a complex representation of the carrier (at zero frequency 
now) and sidebands within fj{2r) of the carrier. Because frames are processed independently, it is all 
right to lose A D data between frames while carrying out further processing on completed z-frames. 



Fig. 2. Frequency response of the FIR filter for lowpass decimation. 


B. Signal Spectrum 

The signal spectrum is obtained as a two-sided spectrum of the z-frame. First, the z-frame is scaled to 
unit energy. Most of the energy is in the carrier, which is now at dc (zero frequency). To prevent this dc 
energy from leaking into the rest of the spectrum, we get rid of most of it by removing a linear fit from 
the frame^ We call this kind of preconditioning operation a calibration. The specific example used here 
canbedefined on a general array y[0], ■ ■ ■ ,y[N - 1] as follows: Let Af be an integer approximately equal 
to iV/6. Compute the centroid points 


l M ~ 1 j /v-i 

(to,yo) = ^ ( n >yM), (ti,yi) = — V' (n,y[n]) 

"=» M n=N — M 


and pass a straight line co + cm through them. The calibrated array is given by t/ 0 fn] = ufnl - cn - cm 
If y itself is a straight line, then y 0 = 0. 1 J W ‘ 00 l 

The choice of this particular operation (especially the N/ 6) for spectral preconditioning is admittedly 
seat-ot-the-pants engineering. Perhaps removing a conventional least-squares fit would do as well. To 
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deal with time series modeled by processes that are possibly nonstationary but do have stationary first 
or second increments, it is desirable to subtract some linear fit, not just the mean. This makes all fram 
statistically identical, so that the average of the spectral estimates of J disjoint frames converges as 
j _ ► oo, just as in the theory of stationary-process spectral estimates. 

The spectrum and frequency arrays are now given in terms of the calibrated array z 0 by 


S z [k] = 


f,N z f |tf r (/[*])l 


N.f-l 


n=0 


( - i27mk\ 
) ] 2o[n]u 0 [n;iV z f,4]exp I — — J 


m = 


rN fft 


where k = -N« t /2 + !,•••, N fft /2. The squared magnitude of 


nh-1 / 

H r {f) = Ys M»]ex P ( 

n=0 ' 


— i2Trnf\ 

fs ) 


is used for equalizing the spectrum against the lowpass decimation filter. As b ^ a 
is labeled dBc/Hz. Points corresponding to frequencies with absolute value below /»/ ( ** ) 

95 percent of the Nyquist frequency 0.5/, /r are not displayed. The low cutoff hides doubtful values 
dc P the high cutoff hides a 3-dB rise at the Nyquist frequency caused by the combination of lowp 
decimation, noise folding at the Nyquist frequency, and equalization. The .^ er ch °°^ S ° W many ° 
these frame spectra are averaged into a run spectrum. The resolution bandwidth is given by 


Wnb = 


r (X>oM) 2 


(4) 


C. Amplitude and Phase 

Extraction of amplitude and phase residuals starts with a rectangular-to-polar operation on the 
z-frame. The result is a complex “amplitude-phase” frame ap[n], n = 0 - 1, whose^reai pa 

is the amplitude of z\n\ and whose imaginary part 8\n] is the phase of z[n] wrapped in [ , ] • 

LpUtute are repjed by their fractional deviations from the mean. The phases are unwrapped mto 
phase deviations 0[n] (replacing 6[n\ in the ap array) by the following algorithm: 


0[O] = 0, <f>[n\ = 4>[n - 1] + mods (0[n] - 6[n - 1],2tt), 


n = N zi - 1 


The symmetric residue function mods is defined by 


mods(x, a) = x -a nint (x/a) W 

The correctness of this algorithm requires only that |A#i]| < tt. The mods function also plays the central 
role in the unwrapping algorithm described in Appendix C. 
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The amplitude or phase residuals can be displayed as time series for the frame. Often, the phase 
residuals are dominated by a ramp (a frequency offset), so that it is desirable to subtract a linear fit 
to reveal the random fluctuations. This can be done with the calibration operation described above in 
connection with spectral preconditioning. 

For better or worse, amplitude and phase spectra are computed together by a single complex FFT 
instead of two real FFTs. The real and imaginary parts of the ap[n] array are calibrated as above, tapered 
by uq [n; N zf , 4], and zero-padded to N m elements. Let AP[fc], k = 0, • ■ • , Nff t - 1 be the complex Fourier 
transform of the resulting array. The transforms of the real amplitude and phase frames are given by 

A[fc] = | (AP[fc) + AP[AT fft - *]*) , «[fc] = i (AP[*J - AP[7V fft - It]*) 

for k = 0, ■ • ■ where A P [ Algt ] is defined to be AP[0]. The one-sided amplitude and phase spectra 

for the frame are given by 


S a [k] = 


fsN z{ |i/p(/[*])|- 




fsN z{ \H r (f[k})\' 


|$[ fc ]| 2 


(6) 


with frequency array f[k] = (f s /(rN m ))k. We apply the same low- and high-frequency cutoffs as we did 
with the medium-band signal spectrum. The absence of a factor of 2 in the scaling factor of Eq. (6) [see 
Eq. (1)] gives a single-sideband presentation of the spectra, so that they can be labeled dBc/Hz when 
converted to dB. If a factor of 2 were present in the numerators, the unit for S# would have to be rad 2 /Hz. 
As before, a number of frame spectra can be averaged into a run spectrum. The resolution bandwidth is 
given by Eq. (4). 


V. Narrow-Band Processing 

This processing mode also assumes that the signal consists of a carrier with weak sidebands. Its 
purpose is to achieve an arbitrarily large reduction in data rate, limited only by the user’s patience. The 
stream of A— D data is reduced to a sequence of average amplitude and phase residuals, the averaging 
time being chosen by the user. The phase residuals from two channels can be combined into a differential 
phase. These streams of band-reduced data can be processed into time series, spectra, or Allan deviations 
(phase or differential phase only). 

A. Amplitude and Phase Extraction 

The stream of A-D data is divided into batches of size IV xb , which must be adjacent for the entire run. 
There is a minimum and maximum batch size (now 200 and 8192). A frame consists of n bf batches, or 
AI x f = n,j,f jV x i, A-D data, where n b f can be any positive integer. Each batch is reduced to one sample of 
average amplitude and phase, and n b f batch samples are averaged to produce a frame sample. The user 
has to choose .V xb and n b f (with the bounds on A' xb enforced by the user interface) to achieve the desired 
reduced sample rate f s /N x f. Unless there are phase tracking problems (see below), the results for a fixed 
frame size should depend little on the number of batches per frame. 

Let us represent the digitized signal by 

x(t ) = A(t) cos <f>(t) 

where $(f) is the total phase, which one can think of as wt + 6 + </(<), where 4>{t) is a phase residual. 
The point is that <£(f) is an intrinsic part of the signal (except for an unknowable additive constant 
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2jrno), while w and <f>(t) trade off with each other. We assume that A(t) and $(<) satisfy two imprecisely 
given conditions, called here the assumptions of small local variations: (1) Over a batch, the fractional 
variations of A(t) from its mean are much less than 1, and (2) over at least two batches, the total phase 
differs by much less than one radian from a first-degree polynomial fit (constant phase offset plus constant 
frequency). Over longer time spans, the phase might deviate from a straight-line fit by many radians. 

Let the batches of a run be indexed by k, k = 0, 1, • - •- Batch k starts at time t k = kN x b/ fa- F° r 
the moment, let t run over the sequence of times t k + n//,,n = 0, • • - , N x \, — 1, in batch k. The Pony 
computation (Parts 1 and 2) of Appendix A is used to estimate the local frequency, amplitude, and phase 
of the batch. It gives o k (radians per cycle), A k , and 8 k such that 

x(t) « A* cos (6 k f s (t-t k )+6 k ) ( 7 ) 

(The sign of d k is taken to be the same as the polarity s p „|.) Write w k = o k f,. With the assumptions of 
small local variations, it turns out that, to first order in these variations, 


A k « A k (8) 

ip k := wjfcfo - t k ) + 9 k « <1* (mod 2 tt) (9) 

where t k , A k , and $ k are the averages of t, A{t), and $(t) over batch k. It is important to note that the 
approximation [Eq. (9)] of ip k to (mod 2?r) is better than the approximation of the phase on the right 
side of Eq. (7) to $(t) because the errors in d k and 8 k tend to compensate each other in just the right 
way. 

The average amplitude residual for batch k is computed by a k = A k /Ao - 1. The computation of 
phase residuals is more delicate. According to Eq. (9), V'fc. to first order, is the average total phase of 
the signal in batch k, modulo 27 r. There are two problems. First, there is the 2n ambiguity. Second, 
we would like to have a phase residual instead of the large total phase. Let us use the initial measured 
frequency ljq and phase Bq to calibrate the total phase to a phase residual 


<j){t) = $(t) — LJo(t — to) — @0 


( 10 ) 


where f now runs over all time beyond the starting time t 0 of the run. Note that 4>{t) depends on the 
calibration parameters i2>o t 8o> so it is not intrinsic. Its average over batch k is 


4>k = ~ d)Q (tfc — to) — 


( 11 ) 


These are the batch phase residuals that we would like to compute. From Eq. (9) it follows that, to first 
order, 


4> k « 1 p k - u>o (tfc - to) -do (mod 2 t r) 
0o ~ 0 (mod 27 t) 


( 12 ) 


To a good approximation, then, we know the 4> k , modulo 27r. Because of the assumption of small local 
variations, we also can predict, with an error < tt, how many radians the average total phase advances 
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from one batch to the next, given its previous behavior. With this information, and with the measured 
assumed to be 0, the 27 t ambiguity can be removed sequentially from all the (j> k by means of a second-order 
unwrapping algorithm given in Appendix C. It is the same algorithm, with different parameters, that is 
used for unwrapping the picket fence time-interval measurements that capture the 1-Hz zero crossings. 

The algorithm also produces a sequence of prediction errors z k that satisfies \z k \ < 7 r. It measures how 
much the current phase differs from what we think it should be, based on the behavior of the previous 
batches. If any \z k \ exceeds a certain threshold, now set at tt/2 , a caution is issued to the user. Perhaps 
the frequency is changing so fast that the assumption of small variations fails for the batch length N xb . 
In effect, the analyzer may be losing phase lock, like a phase-locked loop whose bandwidth is too small. 
If this happens, the user can try decreasing N xb . As mentioned above, the amplitude and phase residual 
averages for a frame are obtained simply by averaging n b f batch values. Thus, if the user has to decrease 
iV xb to keep the analyzer in lock, he can maintain his chosen averaging time by increasing n b f. 

B. Differential Phase 

By differential phase we mean some method of subtracting the phases of two channels that are being 
sampled simultaneously at the same rate. There are two flavors of differential phase processing. In S-S 
or X-X differential phase, it is assumed that both channels (1 and 2) originate at the same RF band and 
are downconverted to the same frequency. In this case, the total phases should not be too far apart, and 
so it makes sense to compute the batch averages 


6$ k = $*(1) - $*(2) - 27rn 0 


where (1) and (2) identify the two channels and no is the integer that makes —tt < ^^0 5: ?r- Applying 
Eq. (11) to both channels, we obtain 


6$k - 4>k( 1) - <M2) + (£o(l) - £o(2)) (t k - f 0 ) + <? 0 (1) - 0 O (2) - 27rn 0 (13) 

which gives the intrinsic quantity 6$ k in terms of measured quantities. 

The original design of the analyzer included a sample-and-hold unit so that channels 1 and 2 could be 
sampled simultaneously. This is no longer the case and, hence, the channel samples have to be interleaved 
at total rate 2 f s through the A-D converter: (1), (2), (1), (2), • • • , where a channel 1 sample is paired with 
the following channel 2 sample. To deal with this situation, we use current batch frequency estimates to 
adjust the total phases of the two channels as if they were sampled halfway between the channel 1 sample 
time and the channel 2 sample time. The phase advance of channel 1 over a delay 1/(4 f s ) is estimated 
35 ^/vidCl)/^/^), where / v id(l) } the current estimate of the analog carrier frequency of channel 1, is 
computed by / v id(l) — /s(^base(l) 4- d k (l)/(2n)). A similar correction of opposite sign is applied to the 
channel 2 total phase. Consequently, a correction 

~2 (^base(l) 4- n base (2)) + — (6fc(l) 4- 0/c(2)) 


has to be added to <5$*. 

In S-X differential phase, channel 1 is downconverted from 2295 MHz (S-band), channel 2 from 
8415 MHz (X-band), or the reverse, and we are required to produce some version of 


S band phase - — (X band phase) 
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In a preliminary design, the analyzer simply computed the nonintrinsic quantity 4> k ( 1) - (3/ll)0fc(2), 
which depends on the initial measured frequencies u>o(0> 1 — 1» • ■ * » 2, and which, if a linear fit is not 
removed , has a random ramp component that depends on these measured frequencies. The current 
design uses a more objective method in which the measured frequencies are replaced by a priori known 
design frequencies w 0 (») = / 3 o 0 (i). These are computed from the user-provided analog offset frequencies 
/ofst(*) by w 0 (») = 2tt (/ ofst (t) - / 3 n base (i)). One can then produce phase residuals <t> k {i) = 4> k (i) + (wo(0 
- w 0 (»))( t fc “ to ) that start at zero but show ramps if the actual channel frequencies differ from the design 
frequencies. S-X differential phase is now just <j> k (l) -(3/ll)<A fe (2), which shows a ramp if the frequencies 
of the S- and X-channels are not related in exactly the right way. In contrast with the S-S or X-X 
situations, the first sample of this differential phase is zero; we are calibrating for frequency only and not 
attempting to measure the absolute synchronization of the two channels. 

As with amplitude and phase, the batch averages of differential phase are combined into frame averages. 


C. Time Series 

The stream of narrow-band samples (frame average amplitude residuals, phase residuals, or differential 
phases) can be collected into a buffer and plotted against time. In the present software, we use a buffer 
management scheme that automatically subsamples the buffer by a factor of 2 when it fills up, crunches it 
to half its size, and begins to accept data at half the previous rate. At any time during the run, the buffer 
contains a record of the entire data stream, subsampled by some power of 2. Because phase residuals 
and differential phases are likely to be dominated by a straight line, we normally apply the calibration 
operation described in Section V.B before plotting them so that random fluctuations can be seen. 


D. Spectrum 

Any of the streams of narrow- band samples can be subjected to the same spectral estimation process. 
Because it takes longer to collect the data arrays, there is incentive to use the narrow-band data more 
efficiently than the medium-band data. In compensation, there is more processor time available per 
A-D sample for expensive postprocessing. We use an unweighted Thomson multitaper spectral estimator 
[10,7 (Chapter 7)] with orthogonal data tapers (trig prolates) computed by the author [5] (Appendix B). 
The user chooses a FFT size a power of 2. At the start of the test, we compute an array of K 
orthogonal data tapers u fc [n; N m , tu], n = 0, • • • , N fft - 1, k = 0, • ■ • , K - i. The value of K depends on w 
and on the sidelobe level we wish to tolerate in the frequency responses of the u k . In the present design, 
w = 4, K = 4. An array of samples i[0], • • • , x[A fft - 1], called a “narrow-band frame” (nbframe), is 
preconditioned by the calibration operation of Section V.B. Then K distinct “eigenspectra” So, • • • , Stf-i 
are computed by applying the tapers and a real FFT, giving 


S k {m] 


N* f 

Afft/a 


Nm - 1 

V x[n}u k [n-,N f f t ,w]exp(-i‘2nnm/N m ) 

71 = 0 


(14) 


with frequency array f[m\ = (f s /{N x f A/fr t ))m, m — 0, • ■ • , Nfft/2. The spectrum of the nbframe is 
computed by averaging the eigenspectra: 


K - i 

S[m] = — 

fc= o 


(15) 


and the overall run spectrum is computed by accumulating and averaging all the nbframe spectra. One 
advantage of this method is that, over smooth regions of the true spectrum, the variance of S[m] is about 
K times smaller than the variance of each Sk[m\. With a single-taper method, variance could be reduced 
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by using shorter nbframes 
the resolution bandwidth. 


or averaging the spectrum over frequency. Either of these methods increases 


To prepare the spectrum for display, we cut off frequencies below 
conversion to dBc/Hz. The resolution bandwidth W nh is given by 


(fs/(N xf N m ))w and do the usual 


1 

^nb 


K-l 


1 “ J 


Jt= 0 


1 

W nb , k 


where 


w nhik = 


fsNfh 

^xf(E^M) 2 


is the resolution bandwidth of S k . Although it is not apparent, IF nb is proportional to l/jV fft ; 
use Nff t to trade off resolution against run length. 


one can 


The user should be aware that the spectral window of this method is not bell shaped but approximately 
rectangular with ripples across the top. If the spectrum has a bright line whose width is of the order of 
° ne T , or less ’ the ‘mage of the line may appear to have four small peaks at the top. These are 
artifacts^ tfie met hod and do not indicate a splitting of the line. (See Appendix B and Fig. 3.) 

In the current version of narrow-band processing, we have achieved bandwidth reduction by unweighted 
averaging: The batch samples of amplitude and phase are, to first order, unweighted averages of these 
quantities, and frame samples are unweighted averages of batch samples. Consequently, a calculated 



FFT BIN 


Fig. 3. Shape of a bright line for narrow- 
band spectrum. 
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spectrum for frequencies between 0 and the Nyquist frequency f,/(2N*) is not, strictly speaking, an 
estimate of the spectrum of the quantity in that frequency range but rather an estimate of the spec 
of the averages of the quantity over an averaging time r 0 = N xf / f a , sampled at rate 1/t 0 . Th p 
suffers from both aliasing and distortion. The Pony method of extracting batch samples of amplitude 
and phase leads inherently to this situation for frames consisting of one batchy The mam decision was 
how to deal with further bandwidth reduction: whether to use a lowpass decimation filter, a bank of such 
filters or simply to extend the situation with unweighted averaging. The advantages of the chosen esign 
ml simplicity consistency, and flexibility in the choice of decimation factor (frame length), which can be 
large enough to exhaust the patience of any user. 

E. Allan Deviation 

The stability analyzer can compute the Allan deviation of frame samples of phase or differential phase 
for an array of averaging times r that are powers of 2 times the frame duration r 0 . It was required 
remove an estimate of linear frequency drift from the results. For a drift estimator we use the simp 
three-ooint estimator suggested by Weiss [11]. Although the basic method is covered in [2] and [3j, we 
rui^Throug^the'compu^ions for a particular value of r = nr 0 . Let the stream of phase samples be 
^o, </>!, -• 8 At a given point in the run, we have accumulated sums of the first and second powers o 
second differences of 0, with stride n, namely, 


m-f 1 

S P =E( A n<M P , P = 1 i 2 

where m > 4. (The author realizes that the sum for p = 1 telescopes.) We have also collected a 
subsampled version 


00 , 0d 5 02d>‘ * * 


of the whole run so far by the same buffer mechanism used for time series above. The calculations proceed 
as follows; 


n c = 


/ 

2 


d 


D c = 02n c - 20n c + 00 (unsealed drift estimate) 


V = 


52 

m 



(sample variance) 


V = V + 



(drift correction) 


v = ( m _ i) (q. 8776 + 0.0643e (1/2)(m 4) ) 


(degrees of freedom) 


V ± = V 
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±/ x 

v T ~ V2 2l Tf re{ T 


(Allan deviation with error bar) 


°v(t) = 


W 


V2 2nf Tef r' 


The formula for degrees of freedom is an empirical formula fitted to the author’s numerical results for the 
random- walk-of-frequency model of phase deviations (/“ 4 noise). The error bars, which are really the 
square roots of one-sigma error bars for <7y(r), should be conservative for f Li noise, (3 > —4. 


VI. Zero-Crossing Processing 

To capture the up-crossing times of the 1-Hz square wave, a preliminary measurement of the nominal 
period p of the square wave is taken with the interval timer, which is then set to measure the time 
intervals between each subsequent up-crossing and the next pulse of a 10-Hz train of reference pulses, 
the “picket fence.” These readings are unwrapped into a sequence of time residuals, as described in 
[4]. The algorithm, which is really the same as the one used for unwrapping the narrow- band phase 
deviations (Appendix C), need not be reproduced here. The time deviations produced by this algorithm 
are multiplied by the scale factor 


2?r / ref 
fmixP 


to give phase deviations that can be used like the batch averages of phase deviation that come from the 
narrow-band process. For time series and Allan deviation, we allow only one batch per frame, as the 
1-second period is natural for the user. For spectrum, an arbitrary number of batches per frame is allowed 
so that users can shrink the Nyquist frequency below 0.5 Hz as much as they want. 
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Appendix A 
The Pony Calculation 


This is a batch method for computing the frequency, amplitude, and phase of a sampled sine wave. 
It comes from a method of harmonic analysis called Prony’s method [6 (Chapter 11)], which analyzes 
a waveform into the sum of n sine waves. The calculation we call “Pony” is simply a modification of 

Prony’s method for n = 1. 


I. Part 1 : Frequency 

Let the data array be x[0\, ■■■,x[N- 1]. If x[n) were exactly of form A cos (on + 6), then we would 
have 


x[n + 1] +x[n - 1] = (2 cos o)x[n], n = 1,- -,N - 2 (A- 1 ) 

On the other hand, if x[n] is a noisy cosine wave, then let us estimate coso by projecting the vector 
x [ n + l] + x [n - 1] orthogonally onto the vector x[n}. The computation is 

(1 /2) ( x[0]s[l] + x[N - 2|4jV - 1]) + EnJi 3 *["]*(" + 1] 

c_ E£i a *l »] 2 ~~ 


o = arccos (c) in [0,7r] 

if |d < 1, else o goes to the nearest port in the storm, 0 or n. One may also change the sign of o according 
to polarity considerations. 
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whKe'o = I"™?-, 2 !™' 1 el5e " h «V 8lnele - prec “ ion com P>“ army Of powers - 0 • • • W - 1 

Compute the dyidic powt<' J TudoXT'’ ^ 7 the ““ " = * 

Multiply by J to " “.4 ^ 7 “T" ,he produc,s *° « 8 , «”■ 

steps ge, more J more 'effldem “or'a USS'^T *“ ** "" S ““™ 


II. Part 2: Amplitude and Phase 


Having estimated the frequency, we 
and solve the least-squares problem 


use it to estimate amplitude and phase. Let a = A cos 0, b 


A sin 9 , 


x[n] % a cos on - b sin on, n = 0, • • • , TV - 1 

JfoZ* cients of the normal matrix can easily be expressed in closed *«». 


w-l 


n=0 


-i[' 


cc = ^ (A^ + cos (o(JV- 1))2? oAr 


sin o 


yv-i 


z[n] cos on, x[n] si 


sin on 


n=0 


ss = - 
2 


AT - cos (o(N - 1))^] 


sin 


(A-2) 


cs= - sin (o(N - 1)) 


sin oN 
sin o 


D = 


cc • ss - cs 


a = ss ’ 4- cs • x s b __ cs • x c 4- cc • x fl 


£> 




-4 — \/a 2 + b 2 , 9 = angle(a + ift) 

Most of the work is in the in-phase and quadrature mixing; ooeratinn fPn (A o\i r,* u l 

« n whose generation is described in Part 1. g ° Peratl ° n ^ whlch uses the array 

The calculation given here can be regarded as an improvement on the approximations 


r, ~ 2 . 2 

~ N Xc ’ b * N x ° 
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Appendix B 
Windows and Filter 


in place of the notation Uk{n; N,w/N,w] in [5]. 

Figure 1 shows the frequency responses (spectral 

tion. The Qos curve applies to full-band spectrum, 04 e igenspectra, E q. (14), that are averaged 

spectra. Note that fi 4 is the ^® rag ® 0 f a spectral estimate is’ the convolution of the true 

^ appear £ £ s^ecTral 

SmaTint IsTp^oZZLT^tZ ripples at the top will not be so prominent on a typical 
dB scale. 

The jV-point FIR lowpass filter used in medium-band processing before decimation by r is b 
coiwentiona^way from the trig prolate window u 0 [n; N,w\. The formula for it is 


= uo[n; N,w ] sine ^27r/h (n - ^ > 


n = 0, • • • , N — 1 


normalized so that 5Z ^r[ n ] ^ where 


w 


-4 N = 16 r, fh = 


0.4 


sinx 


sme x — 


r fv.?c filtpr for r = 2 The response is essentially the same for 
Figure 2 shows the frequency Fig . 2 . Only one table is needed to represent the 

all r if frequency is scaled according to tne x axis uirg j 

ft^uency response fo, the purpose of dualising the medium-band spectra. 
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Appendix C 

Phase Unwrapping Algorithm 


This algorithm produces the narrow-band phase residuals 4> k from the carrier frequency and phase 
estimates extracted from each batch by the Pony calculation. It is assumed that the batches all have 
length iV x b and are adjacent. Recall the definition, Eq. (5), of the mods function. Let the damping 
constant A be a number between 0 and 1. In the following algorithm, ip k is related to ip k of the main text 
by ip' k = t/'fc - o 0 (N xh - l)/2. 


ip' 0 = 0 O 

zo = 0, 4> 0 = 0, q 0 = 0 
For k — 1, 2, • 

Obtain the batch frequency and phase d k , 0 k . 

= (6fc - 6 0 )(iVxb - l)/2 + 0 k ! V'fc is total phase mod 27 t. 
z k = mods {i)' k - xp k _ 1 - 6 0 N xb - q k ~ i, 27 t) ! prediction error. 

If \z k \ > 7t/ 2 (say), then issue caution “losing lock” to user. 

= <j> k _i 4- q k „i + z k ! output phase residual. 
q k = q k - 1 + Xz k ! low pass-filtered A(p k . 


Next k. 

Note that q k , z k satisfy 


Qk = ( 1 ~ A)g fc _i + AA0 fc , z k = - qk- 1 


This says that g* is a lowpass-filtered version of A0*, and z k is a prediction error for A <j> k . The basis of 
the algorithm is (1) the assumption that \z k \ < n and (2) the knowledge of z k modulo 27t, namely, 


z k = A$* - £ 0 A t k - q k -i = A^fc - do^xb - Qk-i (mod 2 tt) 


Any value for A in [0, 1] is meaningful. If 0 < A < 1, then, in effect, a weighted average of previous phase 
advances, with weights (1 — A) n , is used to judge what the current phase advance should be. In the script 
files that drive the software, A has been set to 1/10. This provides some stability against large errors 
while maintaining the ability of the algorithm to follow frequency drifts. 
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Errata 


In “Adaptive Line Enhancers for Fast Acquisition” by H.-G. Yeh and T. M. Nguyen, which appeared in 
The Telecommunications and Data Acquisition Progress Report 42-119 , July September 1994 > November 
15, 1994, the plot in Fig. 14 was incorrectly situated. The correct figure is provided below. 
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Fig. 14. Magnitude of the input data to the ALE, ALEDF, AND ALECA. 
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